function [weights, ranks, varargout] = bestRanks(indicators, scale, varargin)
%BESTRANKS - Rank-optimal weighting of indicators.
%bestRanks calculates rank-optimal weightings of Q indicators for C
%countries yielded by first or second order optimization.
%For further information see 'Rank-optimal weighting or "How to be best in
%the OECD`s Better Life Index?"' (J. Lorenz, C. Brauer, D. Lorenz, 2015).
%
% Syntax:  see Example
%
% Inputs:
%    indicators - C x Q matrix
%    scale - 'int' (integer) or 'cont' (continuous)
%    criterion (optional) - 'distance' or 'score' (second order criterion)
%
% Outputs:
%    weights - C x Q matrix
%    ranks - m vector
%    distances/scores (optional) - m vector
%
% Example: 
%    indicators = randn(20, 10);
%    [weights, ranks] = bestRanks(indicators, 'int');
%    [weights, ranks, distances] = bestRanks(indicators, 'int', 'distance');
%    [weights, ranks, scores] = bestRanks(indicators, 'int', 'score');
%
% Other m-files required: none
% Subfunctions: writeCIP,  writeMAT
% MAT-files required: none
% SCIP!!!
%
% See also: -

% Author: Christoph Brauer
% TU Braunschweig
% email: ch.brauer@tu-bs.de
% Website: https://www.tu-braunschweig.de/iaa/personal/brauer
% July 2015; Last revision: 22-July-2015


    %% validate input arguments
    p = inputParser;
    p.addRequired('scoreMatrix', @isnumeric);
    p.addRequired('scale', @ischar);
    p.addOptional('criterion', '', @ischar);
    p.FunctionName = 'bestRanks';
    p.parse(indicators, scale, varargin{:});
    criterion = p.Results.criterion;
    assert(strcmp(scale, 'int') || strcmp(scale, 'cont'), 'The value of ''scale'' is invalid. It must be ''int'' or ''cont''.');
    assert(strcmp(criterion, 'distance') || strcmp(criterion, 'score') || strcmp(criterion, ''), 'The value of ''criterion'' is invalid. It must  be ''distance'' or ''score''.');

    %% first order step
    
    % determine number of countries and indicators
    [C, Q] = size(indicators);
    
    % write first order problem for each country in CIP format
    writeCIP(indicators, C, Q, scale);
    
    % solve first order problem for each country using SCIP
    for k = 1:C
        system(sprintf('scip -c "read country_%02.0f.cip optimize write solution country_%02.0f.sol quit"', k, k)); 
    end
    
    % write results for first order problem in MAT format
    writeMAT(C, Q);
    
    %% second order step
    if nargin == 3
        
        % write second order problem for each country in CIP format (first
        % order CIP files are overwritten)
        writeCIP(indicators, C, Q, scale, criterion);

        % solve second order problem for each country using SCIP (first
        % order SCIP output files are overwritten)
        for k = 1:C
            system(sprintf('scip -c "read country_%02.0f.cip optimize write solution country_%02.0f.sol quit"', k, k)); 
        end

        % write results for second order problem in matlab format (first order
        % MAT files are overwritten)
        writeMAT(C, Q, criterion);
    end
    
    %% write output
    weights = zeros(C, Q);
    ranks = zeros(C, 1);
    
    % additional output in case of second order problem
    if nargin == 3
        
        % distance maximization
        if strcmp(varargin{1}, 'distance')
            distances = zeros(C, 1);
            
        % score maximization
        else
            scores = zeros(C, 1);
        end
    end
    
    % add results for each country
    for k = 1:C
        
        % load results of country k
        load(sprintf('country_%02.0f.mat', k));
        
        % add results of country k to output
        weights(k, :) = w;
        ranks(k) = bestRank;
        
        % additional output in case of second order problem
        if nargin == 3
            
            % distance maximization
            if strcmp(varargin{1}, 'distance')
                distances(k) = bestDistance;
                
            % score maximization
            else
                scores(k) = bestScore;
            end
        end
    end
    
    % set variable output in case of second order problem
    if nargin == 3
        
        % distance maximization
        if strcmp(varargin{1}, 'distance')
            varargout{1} = distances;
            
        % score maximization
        else
            varargout{1} = scores;
        end
    end
    
    %% remove temporary files
    system('rm *.cip');
    system('rm *.mat');
    system('rm *.sol');
    
end

%% writeCIP
function [] = writeCIP(indicators, C, Q, varargin)

    % first order problem
    if nargin == 4
        
        % integer weights
        if strcmp(varargin{1}, 'int')
            
            % write one CIP file for each country
            for k = 1:C
                
                % calculate matrix for problem constraints
                indicatorDiff = ones(size(indicators, 1) - 1, 1) * indicators(k, :) - [indicators(1:k-1, :); indicators(k+1:end, :)];
                
                % write file
                fStr = sprintf('country_%02.0f.cip', k);
                fId = fopen(fStr, 'w');
                fprintf(fId, 'STATISTICS\n');
                fprintf(fId, 'Problem name\t: %s\n', fStr);
                fprintf(fId, 'Variables\t\t: 81 (%i binary, %i integer, 0 implicit integer, %i continuous)\n', C - 1, Q, C - 1);
                fprintf(fId, 'Constraints\t\t: 71 initial, 71 maximal\n');
                fprintf(fId, 'OBJECTIVE\n');
                fprintf(fId, 'Sense\t\t\t: maximize\n');
                fprintf(fId, 'VARIABLES\n');
                for i = 1 : C - 1
                   fprintf(fId, '\t[binary] <z%i>: obj=1, original bounds=[0,1]\n', i);
                end
                for i = 1 : Q
                   fprintf(fId, '\t[integer] <w%i>: obj=0, original bounds=[0,5]\n', i);
                end
                for i = 1 : C - 1
                   fprintf(fId, '\t[continuous] <s%i>: obj=0, original bounds=[-inf,+inf]\n', i);
                end
                fprintf(fId, 'CONSTRAINTS\n');
                fprintf(fId, '\t[linear] <c0>:');
                for i = 1:Q - 1
                   fprintf(fId, ' <w%i> +', i);
                end
                fprintf(fId, ' <w%i> >= 1;\n', Q);
                for i = 1:C - 1
                   fprintf(fId, '\t[linear] <c%i>: <s%i>', i, i);
                   for j = 1:Q
                      fprintf(fId, ' %+f<w%i>', -indicatorDiff(i, j), j); 
                   end
                   fprintf(fId, ' == 0;\n');
                end
                for i = 1:C - 1
                    fprintf(fId, '\t[nonlinear] <qc%i>: (<s%i> * <z%i>) >= 0;\n', i, i, i);
                end
                fprintf(fId, 'END\n');
                fclose(fId);
            end
            
        % continuous weights
        else
            
            % write one CIP file for each country
            for k = 1:C

                % calculate matrix for problem constraints
                indicatorDiff = ones(size(indicators, 1) - 1, 1) * indicators(k, :) - [indicators(1:k-1, :); indicators(k+1:end, :)];

                % write file
                fStr = sprintf('country_%02.0f.cip', k);
                fId = fopen(fStr, 'w');
                fprintf(fId, 'STATISTICS\n');
                fprintf(fId, 'Problem name\t: %s\n', fStr);
                fprintf(fId, 'Variables\t\t: 81 (%i binary, 0 integer, 0 implicit integer, 46 continuous)\n', C - 1);
                fprintf(fId, 'Constraints\t\t: 71 initial, 71 maximal\n');
                fprintf(fId, 'OBJECTIVE\n');
                fprintf(fId, 'Sense\t\t\t: maximize\n');
                fprintf(fId, 'VARIABLES\n');
                for i = 1 : C - 1
                   fprintf(fId, '\t[binary] <z%i>: obj=1, original bounds=[0,1]\n', i);
                end
                for i = 1 : Q
                   fprintf(fId, '\t[continuous] <w%i>: obj=0, original bounds=[0,+inf]\n', i);
                end
                for i = 1 : C - 1
                   fprintf(fId, '\t[continuous] <s%i>: obj=0, original bounds=[-inf,+inf]\n', i);
                end
                fprintf(fId, 'CONSTRAINTS\n');
                fprintf(fId, '\t[linear] <c0>:');
                for i = 1:Q - 1
                   fprintf(fId, ' <w%i> +', i);
                end
                fprintf(fId, ' <w%i> == 10;\n', Q);
                for i = 1:C - 1
                   fprintf(fId, '\t[linear] <c%i>: <s%i>', i, i);
                   for j = 1:Q
                      fprintf(fId, ' %+f<w%i>', -indicatorDiff(i, j), j); 
                   end
                   fprintf(fId, ' == 0;\n');
                end
                for i = 1:C - 1
                    fprintf(fId, '\t[nonlinear] <qc%i>: (<s%i> * <z%i>) >= 0;\n', i, i, i);
                end
                fprintf(fId, 'END\n');
                fclose(fId);
            end
        end
        
    % second order problem
    else
        
        % integer weights
        if strcmp(varargin{1}, 'int')
            
            % distance maximization
            if strcmp(varargin{2}, 'distance')

                % write one CIP file for each country
                for k = 1:C

                    % calculate matrix for problem constraints
                    indicatorDiff = ones(size(indicators, 1) - 1, 1) * indicators(k, :) - [indicators(1:k-1, :); indicators(k+1:end, :)];

                    % load results of first optimization step
                    load(sprintf('country_%02.0f.mat', k));

                    % write file
                    fStr = sprintf('country_%02.0f.cip', k);
                    fId = fopen(fStr, 'w');
                    fprintf(fId, 'STATISTICS\n');
                    fprintf(fId, 'Problem name\t: %s\n', fStr);
                    fprintf(fId, 'Variables\t\t: 83 (%i binary, %i integer, 0 implicit integer, %i continuous)\n', C - 1, Q + 1, C);
                    fprintf(fId, 'Constraints\t\t: 73 initial, 73 maximal\n');
                    fprintf(fId, 'OBJECTIVE\n');
                    fprintf(fId, 'Sense\t\t\t: maximize\n');
                    fprintf(fId, 'VARIABLES\n');
                    for i = 1 : C - 1
                       fprintf(fId, '\t[binary] <z%i>: obj=0, original bounds=[0,1]\n', i);
                    end
                    for i = 1 : Q
                       fprintf(fId, '\t[integer] <w%i>: obj=0, original bounds=[0,5]\n', i);
                    end
                    fprintf(fId, '\t[integer] <sumW>: obj=0, original bounds=[0,+inf]\n');
                    for i = 1 : C - 1
                       fprintf(fId, '\t[continuous] <s%i>: obj=0, original bounds=[-inf,+inf]\n', i);
                    end
                    fprintf(fId, '\t[continuous] <d>: obj=1, original bounds=[0,+inf]\n');
                    fprintf(fId, 'CONSTRAINTS\n');
                    fprintf(fId, '\t[linear] <c1>: <w1>');
                    for i = 2:Q
                       fprintf(fId, ' +<w%i>', i); 
                    end
                    fprintf(fId, ' -<sumW> == 0;\n');
                    fprintf(fId, '\t[linear] <c2>: <sumW> >= 1;\n');
                    fprintf(fId, '\t[linear] <c3>: <z1>');
                    for i = 2:C - 1
                        fprintf(fId, ' +<z%i>', i);
                    end
                    fprintf(fId, ' == %i;\n', C - bestRank);
                    for i = 1:C - 1
                       fprintf(fId, '\t[nonlinear] <qca%i>: <s%i>', i, i);
                       for j = 1:Q
                          fprintf(fId, ' %+f<w%i>', -indicatorDiff(i, j), j); 
                       end
                       fprintf(fId, ' +0.1(<d> * <sumW>) == 0;\n');
                    end
                    for i = 1:C - 1
                        fprintf(fId, '\t[nonlinear] <qcb%i>: (<s%i> * <z%i>) >= 0;\n', i, i, i);
                    end
                    fprintf(fId, 'END\n');
                    fclose(fId);
                end
                
            % score maximization
            else
                
                % write one CIP file for each country
                for k = 1:C

                    % calculate matrix for problem constraints
                    indicatorDiff = ones(size(indicators, 1) - 1, 1) * indicators(k, :) - [indicators(1:k-1, :); indicators(k+1:end, :)];

                    % load results of first optimization step
                    load(sprintf('country_%02.0f.mat', k));

                    % write file
                    fStr = sprintf('country_%02.0f.cip', k);
                    fId = fopen(fStr, 'w');
                    fprintf(fId, 'STATISTICS\n');
                    fprintf(fId, 'Problem name\t: %s\n', fStr);
                    fprintf(fId, 'Variables\t\t: 81 (%i binary, %i integer, 0 implicit integer, %i continuous)\n', C - 1, Q + 1, C);
                    fprintf(fId, 'Constraints\t\t: 74 initial, 74 maximal\n');
                    fprintf(fId, 'OBJECTIVE\n');
                    fprintf(fId, 'Sense\t\t\t: maximize\n');
                    fprintf(fId, 'VARIABLES\n');
                    for i = 1 : C - 1
                       fprintf(fId, '\t[binary] <z%i>: obj=0, original bounds=[0,1]\n', i);
                    end
                    for i = 1 : Q
                       fprintf(fId, '\t[integer] <w%i>: obj=0, original bounds=[0,5]\n', i);
                    end
                    fprintf(fId, '\t[integer] <sumW>: obj=0, original bounds=[0,+inf]\n');
                    for i = 1 : C - 1
                       fprintf(fId, '\t[continuous] <s%i>: obj=0, original bounds=[-inf,+inf]\n', i);
                    end
                    fprintf(fId, '\t[continuous] <score>: obj=1, original bounds=[0,+inf]\n');

                    fprintf(fId, 'CONSTRAINTS\n');
                    fprintf(fId, '\t[linear] <c0>:');
                    for i = 1:Q - 1
                       fprintf(fId, ' <w%i> +', i);
                    end
                    fprintf(fId, ' <w%i> >= 1;\n', Q);
                    for i = 1:C - 1
                       fprintf(fId, '\t[linear] <c%i>: <s%i>', i, i);
                       for j = 1:Q
                          fprintf(fId, ' %+f<w%i>', -indicatorDiff(i, j), j); 
                       end
                       fprintf(fId, ' == 0;\n');
                    end
                    for i = 1:C - 1
                        fprintf(fId, '\t[nonlinear] <qc%i>: (<s%i> * <z%i>) >= 0;\n', i, i, i);
                    end
                    fprintf(fId, '\t[linear] <c%i>: <z1>', C);
                    for i = 2:C - 1
                        fprintf(fId, ' + <z%i>', i);
                    end
                    fprintf(fId, ' == %i;\n', C - bestRank);
                    fprintf(fId, '\t[linear] <c37>: <w1>');
                    for i = 2:Q
                       fprintf(fId, ' +<w%i>', i); 
                    end
                    fprintf(fId, ' -<sumW> == 0;\n');
                    fprintf(fId, '\t[nonlinear] <qc%i>:', C);
                    for i = 1:Q
                       fprintf(fId, ' %+f<w%i>', -indicators(k, i), i);
                    end
                    fprintf(fId, ' +0.1(<score> * <sumW>) == 0;\n');
                    fprintf(fId, 'END\n');
                    fclose(fId);
                end
            end
            
        % continuous weights
        else
            
            % distance maximization
            if strcmp(varargin{2}, 'distance')
            
                % write one CIP file for each country
                for k = 1:C

                    % calculate matrix for problem constraints
                    indicatorDiff = ones(size(indicators, 1) - 1, 1) * indicators(k, :) - [indicators(1:k-1, :); indicators(k+1:end, :)];

                    % load results of first optimization step
                    load(sprintf('country_%02.0f.mat', k));

                    % write file
                    fStr = sprintf('country_%02.0f.cip', k);
                    fId = fopen(fStr, 'w');
                    fprintf(fId, 'STATISTICS\n');
                    fprintf(fId, 'Problem name\t: %s\n', fStr);
                    fprintf(fId, 'Variables\t\t: 82 (%i binary, 0 integer, 0 implicit integer, 47 continuous)\n', C - 1);
                    fprintf(fId, 'Constraints\t\t: 72 initial, 72 maximal\n');
                    fprintf(fId, 'OBJECTIVE\n');
                    fprintf(fId, 'Sense\t\t\t: maximize\n');
                    fprintf(fId, 'VARIABLES\n');
                    for i = 1 : C - 1
                       fprintf(fId, '\t[binary] <z%i>: obj=0, original bounds=[0,1]\n', i);
                    end
                    for i = 1 : Q
                       fprintf(fId, '\t[continuous] <w%i>: obj=0, original bounds=[0,+inf]\n', i);
                    end
                    for i = 1 : C - 1
                       fprintf(fId, '\t[continuous] <s%i>: obj=0, original bounds=[-inf,+inf]\n', i);
                    end
                    fprintf(fId, '\t[continuous] <d>: obj=1, original bounds=[0,+inf]\n');
                    fprintf(fId, 'CONSTRAINTS\n');
                    fprintf(fId, '\t[linear] <c0>:');
                    for i = 1:Q - 1
                       fprintf(fId, ' <w%i> +', i);
                    end
                    fprintf(fId, ' <w%i> == 10;\n', Q);
                    for i = 1:C - 1
                       fprintf(fId, '\t[linear] <c%i>: <s%i>', i, i);
                       for j = 1:Q
                          fprintf(fId, ' %+f<w%i>', -indicatorDiff(i, j), j); 
                       end
                       fprintf(fId, ' +<d> == 0;\n');
                    end
                    for i = 1:C - 1
                        fprintf(fId, '\t[nonlinear] <qc%i>: (<s%i> * <z%i>) >= 0;\n', i, i, i);
                    end
                    fprintf(fId, '\t[linear] <c%i>: <z1>', C);
                    for i = 2:C - 1
                        fprintf(fId, ' + <z%i>', i);
                    end
                    fprintf(fId, ' == %i;\n', C - bestRank);
                    fprintf(fId, 'END\n');
                    fclose(fId);
                end
                
            % score maximization
            else
                
                % write one CIP file for each country
                for k = 1:C
                    
                    % calculate constraint matrix
                    indicatorDiff = ones(size(indicators, 1) - 1, 1) * indicators(k, :) - [indicators(1:k-1, :); indicators(k+1:end, :)];

                    % load results of first optimization step
                    load(sprintf('country_%02.0f.mat', k));

                    % write problem into .cip file
                    fStr = sprintf('country_%02.0f.cip', k);
                    fId = fopen(fStr, 'w');
                    fprintf(fId, 'STATISTICS\n');
                    fprintf(fId, 'Problem name\t: %s\n', fStr);
                    fprintf(fId, 'Variables\t\t: 81 (%i binary, 0 integer, 0 implicit integer, 46 continuous)\n', C - 1);
                    fprintf(fId, 'Constraints\t\t: 72 initial, 72 maximal\n');
                    fprintf(fId, 'OBJECTIVE\n');
                    fprintf(fId, 'Sense\t\t\t: maximize\n');
                    fprintf(fId, 'VARIABLES\n');
                    for i = 1 : C - 1
                        fprintf(fId, '\t[binary] <z%i>: obj=0, original bounds=[0,1]\n', i);
                    end
                    for i = 1 : Q
                        fprintf(fId, '\t[continuous] <w%i>: obj=%+f, original bounds=[0,+inf]\n', i, indicators(k, i));
                    end
                    for i = 1 : C - 1
                        fprintf(fId, '\t[continuous] <s%i>: obj=0, original bounds=[-inf,+inf]\n', i);
                    end
                    fprintf(fId, 'CONSTRAINTS\n');
                    fprintf(fId, '\t[linear] <c0>:');
                    for i = 1:Q - 1
                        fprintf(fId, ' <w%i> +', i);
                    end
                    fprintf(fId, ' <w%i> == 10;\n', Q);
                    for i = 1:C - 1
                        fprintf(fId, '\t[linear] <c%i>: <s%i>', i, i);
                        for j = 1:Q
                            fprintf(fId, ' %+f<w%i>', -indicatorDiff(i, j), j); 
                        end
                        fprintf(fId, ' == 0;\n');
                    end
                    for i = 1:C - 1
                        fprintf(fId, '\t[nonlinear] <qc%i>: (<s%i> * <z%i>) >= 0;\n', i, i, i);
                    end
                    fprintf(fId, '\t[linear] <c%i>: <z1>', C);
                    for i = 2:C - 1
                        fprintf(fId, ' + <z%i>', i);
                    end
                    fprintf(fId, ' == %i;\n', C - bestRank);
                    fprintf(fId, 'END\n');
                    fclose(fId);
                end
            end
        end
    end
end

%% writeMAT
function [] = writeMAT(C, Q, varargin)

    % first order problem 
    if nargin == 2
        
        % write one MAT file for each country
        for k = 1 : C

            % open scip output file
            fStr = sprintf('country_%02.0f', k);
            fId = fopen(sprintf('%s.sol', fStr), 'r');

            % parse file
            content = textscan(fId, '%s');
            size = numel(content{1});
            w = zeros(1, Q);
            bestRank = 0; %#ok<*NASGU>
            for i = 1 : size
               if content{1}{i}(1) == 'w'
                   index = str2double(content{1}{i}(2:end));
                   value = str2double(content{1}{i+1});
                   w(index) = value;
               elseif content{1}{i}(1) == 'o'
                   bestRank = C - str2double(content{1}{i+2});
               end
            end
            fclose(fId);

            % save best rank and corresponding weights
            save(sprintf('%s.mat', fStr), 'w', 'bestRank');
        end     
        
    % second order problem
    else
        % distance maximization
        if strcmp(varargin{1}, 'distance')
            
            % write one MAT file for each country
            for k = 1 : C

                % open scip output file
                fStr = sprintf('country_%02.0f', k);
                fId = fopen(sprintf('%s.sol', fStr), 'r');

                % parse file
                content = textscan(fId, '%s');
                size = numel(content{1});
                w = zeros(1, Q);
                bestRank = 0;
                bestDistance = 0;
                for i = 1 : size
                    if content{1}{i}(1) == 'w'
                        index = str2double(content{1}{i}(2:end));
                        value = str2double(content{1}{i+1});
                        w(index) = value;
                    elseif content{1}{i}(1) == 'z'
                        bestRank = bestRank + str2double(content{1}{i+1});
                    elseif content{1}{i}(1) == 'o'
                        bestDistance = str2double(content{1}{i+2});
                    end
                end
                bestRank = C - bestRank;
                fclose(fId);

                % save best rank, distance and corresponding weights
                save(sprintf('%s.mat', fStr), 'w', 'bestRank', 'bestDistance');
            end
         
         % score maximization 
        else
            
            % write one MAT file for each country
            for k = 1 : C

                % open scip output file
                fStr = sprintf('country_%02.0f', k);
                fId = fopen(sprintf('%s.sol', fStr), 'r');

                % parse file
                content = textscan(fId, '%s');
                size = numel(content{1});
                w = zeros(1, Q);
                bestRank = 0;
                bestScore = 0;
                for i = 1 : size
                    if content{1}{i}(1) == 'w'
                        index = str2double(content{1}{i}(2:end));
                        value = str2double(content{1}{i+1});
                        w(index) = value;
                    elseif content{1}{i}(1) == 'z'
                        bestRank = bestRank + str2double(content{1}{i+1});
                    elseif content{1}{i}(1) == 'o'
                        bestScore = str2double(content{1}{i+2});
                    end
                end
                bestRank = C - bestRank;
                fclose(fId);
                
                % save best rank, score and corresponding weights
                save(sprintf('%s.mat', fStr), 'w', 'bestRank', 'bestScore');
            end            
        end
    end
end
